Comprehensive Simulation Results Analysis - All Endpoints

Overview

This document provides a comprehensive analysis of simulation results across all four endpoint types:

• Time-to-Event (TTE): Survival data with Average Hazard Ratio (AHR)
• Binary: Binary outcomes with Odds Ratio (OR)
• Count: Count data with Rate Ratio (RR)
• Continuous: Continuous outcomes with Mean Difference

All analyses are performed consistently across endpoints, with metrics standardized on the log scale (or natural scale for continuous) to enable fair comparisons.

---

Configuration and Data Loading

# --- DEFINE ALL ENDPOINT PARAMETERS ---
ENDPOINT_CONFIGS <- list(
  tte = list(
    folder = "TTE",
    truth_list_name = "simul_parameter",
    truth_overall_name = "true_overall_results",
    truth_subgroup_name = "true_subgroup_ahr",
    truth_metric_exp = "AHR",
    est_metric = "AHR",
    est_col = "estimate_ahr",
    ci_low = "lower_ci_ahr",
    ci_up = "upper_ci_ahr",
    log_truth_col = "true_log_ahr",
    display_name = "Time-to-Event (Survival)",
    metric_scale = "log"
  ),
  binary = list(
    folder = "Binary",
    truth_list_name = "simul_truth",
    truth_overall_name = "truth_overall",
    truth_subgroup_name = "truth_subgroup",
    truth_metric_exp = "exp_value",
    est_metric = "OR",
    est_col = "estimate",
    ci_low = "lower_ci_ahr",
    ci_up = "upper_ci_ahr",
    log_truth_col = "value",
    display_name = "Binary Outcome",
    metric_scale = "log"
  ),
  count = list(
    folder = "Count",
    truth_list_name = "simul_truth",
    truth_overall_name = "truth_overall",
    truth_subgroup_name = "truth_subgroup",
    truth_metric_exp = "exp_value",
    est_metric = "RR",
    est_col = "estimate",
    ci_low = "lower_ci",
    ci_up = "upper_ci",
    log_truth_col = "value",
    display_name = "Count Data",
    metric_scale = "log"
  ),
  continuous = list(
    folder = "Continuous",
    truth_list_name = "simul_truth",
    truth_overall_name = "truth_overall",
    truth_subgroup_name = "truth_subgroup",
    truth_metric_exp = "value",
    est_metric = "Mean Difference",
    est_col = "estimate",
    ci_low = "lower_ci",
    ci_up = "upper_ci",
    log_truth_col = "value",
    display_name = "Continuous Outcome",
    metric_scale = "natural"
  )
)

# --- HELPER FUNCTIONS ---

# Load, validate, and add metadata to RDS files
load_and_validate_rds <- function(file_path, endpoint_id) {
  message(paste("  Loading:", basename(file_path)))
  
  df <- readRDS(file_path)
  
  # VALIDATION: Check for 6000 simulations (6 scenarios * 1000 reps)
  if ("scenario_id" %in% names(df)) {
    n_sims <- df %>% distinct(scenario_id, replication_id) %>% nrow()
  } else if ("scenario_no" %in% names(df)) {
    n_sims <- df %>% distinct(scenario_no, simul_no) %>% nrow()
  } else {
    warning(paste("Could not validate file:", basename(file_path)))
    return(NULL)
  }
  
  if (n_sims != 6000) {
    warning(sprintf("File %s is incomplete! Expected 6000 sims, found %d.", 
                    basename(file_path), n_sims))
  } else {
    message(paste("    ...Validation OK: 6000 simulations found."))
  }
  
  # Add metadata
  df %>% mutate(
    file_name = basename(file_path),
    estimator_name = str_remove(file_name, paste0("^", endpoint_id, "_")) %>%
                     str_remove("\\.rds$")
  )
}

# Standardize truth data
standardize_truth <- function(endpoint_id, ep_params, truth_env) {
  if (endpoint_id == "tte" && exists("simul_parameter", envir = truth_env)) {
    # TTE Truth
    simul_parameter <- get("simul_parameter", envir = truth_env)
    
    real_params_tidy <- simul_parameter[[ep_params$truth_subgroup_name]] %>%
      mutate(scenario_no = as.character(row_number())) %>%
      pivot_longer(
        cols = -scenario_no,
        names_to = "subgroup_real",
        values_to = "truth_exp"
      ) %>%
      mutate(
        join_key = gsub("x_|\\.", "", subgroup_real),
        truth_log = log(truth_exp)
      ) %>%
      dplyr::select(scenario_no, join_key, truth_log, truth_exp)
    
    true_population_results <- simul_parameter[[ep_params$truth_overall_name]] %>%
      mutate(
        scenario_no = as.character(1:6),
        truth_population_exp = !!sym(ep_params$truth_metric_exp),
        truth_population_log = log(truth_population_exp)
      ) %>%
      dplyr::select(scenario_no, truth_population_exp, truth_population_log)
    
  } else if (exists("simul_truth", envir = truth_env)) {
    # Other Endpoints
    simul_truth <- get("simul_truth", envir = truth_env)
    
    real_params_tidy <- simul_truth[[ep_params$truth_subgroup_name]] %>%
      mutate(
        join_key = gsub("x_|\\.|aa|ab|ba|bb", "", subgroup),
        scenario_no = as.character(scenario)
      ) %>%
      rename(truth_value = !!sym(ep_params$log_truth_col)) %>%
      mutate(
        truth_log = if (ep_params$metric_scale == "natural") truth_value else truth_value,
        truth_exp = if (ep_params$metric_scale == "natural") truth_value else exp(truth_value)
      ) %>%
      dplyr::select(scenario_no, join_key, truth_log, truth_exp, truth_value)
    
    true_population_results <- simul_truth[[ep_params$truth_overall_name]] %>%
      mutate(
        scenario_no = as.character(scenario),
        truth_pop_value = !!sym(ep_params$log_truth_col)
      ) %>%
      mutate(
        truth_population_log = if (ep_params$metric_scale == "natural") truth_pop_value else truth_pop_value,
        truth_population_exp = if (ep_params$metric_scale == "natural") truth_pop_value else exp(truth_pop_value)
      ) %>%
      dplyr::select(scenario_no, truth_population_exp, truth_population_log)
  } else {
    stop("Could not find truth data in loaded .RData file.")
  }
  
  list(
    subgroup = real_params_tidy,
    population = true_population_results
  )
}

# Standardize results data
standardize_results <- function(df, endpoint_id, ep_params) {
  estimator_type <- df$estimator_name[1]
  
  if (estimator_type %in% c("population", "subgroup")) {
    # Naive/Population results
    df %>%
      mutate(
        join_key = gsub("S_", "", subgroup) %>% gsub("Overall", "Overall", .),
        scenario_id = as.integer(scenario_no),
        replication_id = as.integer(simul_no),
        estimate_log = estimate,
        ci_lower = lower_ci,
        ci_up = upper_ci,
        estimator_name = estimator
      ) %>%
      dplyr::select(scenario_id, replication_id, estimator_name, join_key, 
                   estimate_log, ci_lower, ci_up)
  } else {
    # Global/OVAT results
    df %>%
      filter(Subgroup != "Overall") %>%
      mutate(
        join_key = gsub("x_", "", Subgroup),
        estimator_name = paste(model_type, prior_name, sep = "_")
      ) %>%
      mutate(join_key = gsub(": ", "", join_key)) %>%
      rename(
        estimate_log = Median,
        ci_lower = CI_Lower,
        ci_up = CI_Upper
      ) %>%
      {
        if (endpoint_id %in% c("tte", "count")) {
          mutate(.,
            estimate_log = log(estimate_log),
            ci_lower = log(ci_lower),
            ci_up = log(ci_up)
          )
        } else {
          .
        }
      } %>%
      dplyr::select(scenario_id, replication_id, estimator_name, join_key, 
                   estimate_log, ci_lower, ci_up)
  }
}

# Load all data for a single endpoint
load_endpoint_data <- function(endpoint_id, ep_params) {
  message(paste("\n=== Loading data for endpoint:", ep_params$display_name, "==="))
  
  # Load results files
  results_path <- file.path(ep_params$folder, "Results")
  all_files <- list.files(results_path, pattern = "\\.rds$", full.names = TRUE)
  
  all_results_list <- lapply(all_files, function(f) {
    load_and_validate_rds(f, endpoint_id)
  })
  
  # Load truth data
  truth_file <- file.path(ep_params$folder, "Scenarios", "truth.RData")
  truth_env <- new.env()
  load(truth_file, envir = truth_env)
  message(paste("  Loaded truth data from:", truth_file))
  
  truth_data <- standardize_truth(endpoint_id, ep_params, truth_env)
  
  # Standardize all results
  standardized_results_list <- lapply(all_results_list, function(df) {
    if (!is.null(df)) {
      standardize_results(df, endpoint_id, ep_params)
    }
  })
  
  # Combine and merge
  all_estimators_df <- bind_rows(standardized_results_list)
  all_estimators_df$scenario_no <- as.character(all_estimators_df$scenario_id)
  # Normalize join keys: remove dots/underscores/spaces and force lower-case
  all_estimators_df$join_key <- all_estimators_df$join_key %>%
    stringr::str_replace_all("[\\._\\s]", "") %>%
    tolower()
  
  merged_df <- left_join(all_estimators_df, truth_data$subgroup, 
                         by = c("scenario_no", "join_key"))
  merged_df_with_pop <- left_join(merged_df, truth_data$population, 
                                  by = "scenario_no")
  
  message(paste("  Data loaded and merged successfully for", ep_params$display_name))
  
  list(
    merged_df = merged_df,
    merged_df_with_pop = merged_df_with_pop,
    endpoint_id = endpoint_id,
    ep_params = ep_params
  )
}

# Calculate performance metrics
calculate_performance <- function(merged_df) {
  merged_df %>%
    group_by(scenario_no, estimator_name) %>%
    summarise(
      rmse = sqrt(mean((estimate_log - truth_log)^2, na.rm = TRUE)),
      mae = mean(abs(estimate_log - truth_log), na.rm = TRUE),
      bias = mean(estimate_log - truth_log, na.rm = TRUE),
      median_ae = median(abs(estimate_log - truth_log), na.rm = TRUE),
      .groups = 'drop'
    )
}

# Calculate subgroup-level performance
calculate_subgroup_performance <- function(merged_df_with_pop, ep_params) {
  threshold <- if (ep_params$metric_scale == "natural") 0.1 else log(1.1)
  
  merged_df_with_pop %>%
    group_by(scenario_no, join_key, estimator_name) %>%
    summarize(
      truth_log = first(truth_log),
      truth_exp = first(truth_exp),
      truth_population_log = first(truth_population_log),
      bias = mean(estimate_log - truth_log, na.rm = TRUE),
      RMSE = sqrt(mean((estimate_log - truth_log)^2, na.rm = TRUE)),
      coverage = mean((truth_log >= ci_lower) & (truth_log <= ci_up), na.rm = TRUE),
      .groups = 'drop'
    ) %>%
    mutate(
      heterogeneous = (abs(truth_log - truth_population_log) > threshold)
    ) %>%
    mutate(coverage = if_else(is.nan(coverage), NA_real_, coverage))
}

# --- LOAD DATA FOR ALL ENDPOINTS ---
all_endpoint_data <- map(names(ENDPOINT_CONFIGS), function(ep_id) {
  load_endpoint_data(ep_id, ENDPOINT_CONFIGS[[ep_id]])
})
names(all_endpoint_data) <- names(ENDPOINT_CONFIGS)

message("\n=== All endpoint data loaded successfully ===\n")

---

Cross-Endpoint Comparison

Overall RMSE Across All Endpoints

This section compares performance metrics across all four endpoint types in a single view.

# Calculate performance for all endpoints
all_performance <- map_dfr(names(all_endpoint_data), function(ep_id) {
  data <- all_endpoint_data[[ep_id]]
  perf <- calculate_performance(data$merged_df)
  perf %>% 
    mutate(endpoint = data$ep_params$display_name) %>%
    filter(estimator_name != 'Global_r2d2_mid_unshrunk_x4')
})

# Scenario labels
scenario_labels <- c(
  `1` = "Positive\n(homogeneous)",
  `2` = "Positive\n(except 1 subgroup)",
  `3` = "Negative\n(heterogeneity)",
  `4` = "Mild\nheterogeneity",
  `5` = "Strong\nheterogeneity",
  `6` = "Misspecified"
)

# Plot RMSE comparison across endpoints
plot_data_all <- all_performance %>%
  mutate(scenario_name = factor(scenario_no, 
                               levels = names(scenario_labels), 
                               labels = scenario_labels))

ggplot(plot_data_all, 
       aes(x = scenario_name, y = rmse, 
           group = estimator_name, color = estimator_name)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  facet_wrap(~endpoint, ncol = 2, scales = "free_y") +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "Overall RMSE Across All Endpoints and Scenarios",
    subtitle = "Performance comparison of estimators across endpoint types",
    y = "RMSE",
    x = "Simulation Scenario",
    color = "Estimator"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    axis.title = element_text(size = 12),
    legend.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Standardized RMSE (Relative to Subgroup Estimator)

# Calculate standardized RMSE for all endpoints
all_standardized <- map_dfr(names(all_endpoint_data), function(ep_id) {
  data <- all_endpoint_data[[ep_id]]
  perf <- calculate_performance(data$merged_df) %>%
    filter(estimator_name != 'Global_r2d2_mid_unshrunk_x4')
  
  # Get baseline (subgroup) RMSE
  baseline <- perf %>%
    filter(estimator_name == "subgroup") %>%
    dplyr::select(scenario_no, rmse_baseline = rmse)
  
  perf %>%
    left_join(baseline, by = "scenario_no") %>%
    mutate(
      standardized_rmse = rmse / rmse_baseline,
      endpoint = data$ep_params$display_name
    )
})

# Plot standardized RMSE
plot_data_std <- all_standardized %>%
  mutate(scenario_name = factor(scenario_no, 
                               levels = names(scenario_labels), 
                               labels = scenario_labels))

ggplot(plot_data_std, 
       aes(x = scenario_name, y = standardized_rmse, 
           group = estimator_name, color = estimator_name)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  geom_hline(yintercept = 1.0, linetype = "dashed", color = "gray30") +
  facet_wrap(~endpoint, ncol = 2) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "Standardized RMSE Across All Endpoints",
    subtitle = "RMSE relative to 'subgroup (no shrinkage)' estimator",
    y = "Standardized RMSE (relative to subgroup estimator)",
    x = "Simulation Scenario",
    color = "Estimator"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    axis.title = element_text(size = 12),
    legend.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

---

Endpoint-Specific Analyses

# Helper functions for tables
display <- function(x, digits = 2) {
  formatC(x, format = "f", digits = digits)
}

getPercentage <- function(x, digits = 0) {
  val <- 100 * x
  ifelse(is.na(val),
    "NA",
    paste(formatC(val, digits = digits, format = "f"), "%", sep = "")
  )
}

# Generate sections for each endpoint
for (ep_id in names(all_endpoint_data)) {
  ep_data <- all_endpoint_data[[ep_id]]
  ep_params <- ep_data$ep_params
  
  cat("\n## ", ep_params$display_name, " (", ep_params$est_metric, ")\n\n")
  
  # Calculate metrics
  performance_results <- calculate_performance(ep_data$merged_df)
  rmse_results <- performance_results %>% 
    dplyr::select(scenario_no, estimator_name, rmse) %>%
    filter(estimator_name != 'Global_r2d2_mid_unshrunk_x4')
  
  # Plot 1: Overall RMSE
  cat("\n### Overall RMSE by Scenario\n\n")
  
  plot_data <- rmse_results %>%
    mutate(scenario_name = factor(scenario_no, 
                                 levels = names(scenario_labels), 
                                 labels = scenario_labels))
  
  p1 <- ggplot(plot_data,
               aes(x = scenario_name, y = rmse, 
                   group = estimator_name, color = estimator_name)) +
    geom_line(size = 1.2) +
    geom_point(size = 3) +
    scale_color_brewer(palette = "Set1") +
    labs(
      title = paste("Overall RMSE -", ep_params$display_name),
      subtitle = paste("Performance of various estimators of", ep_params$est_metric),
      y = "RMSE",
      x = "Simulation Scenario",
      color = "Estimator"
    ) +
    theme_bw() +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      axis.title = element_text(size = 11),
      legend.title = element_text(face = "bold"),
      panel.grid.minor = element_blank()
    )
  
  print(p1)
  cat("\n\n")
  
  # Plot 2: Standardized RMSE
  cat("\n### Standardized RMSE (Relative to Subgroup Estimator)\n\n")
  
  baseline_rmse <- rmse_results %>%
    filter(estimator_name == "subgroup") %>%
    dplyr::select(scenario_no, rmse_baseline = rmse)
  
  plot_data_std <- rmse_results %>%
    left_join(baseline_rmse, by = "scenario_no") %>%
    mutate(
      standardized_rmse = rmse / rmse_baseline,
      scenario_name = factor(scenario_no, 
                           levels = names(scenario_labels), 
                           labels = scenario_labels)
    )
  
  p2 <- ggplot(plot_data_std,
               aes(x = scenario_name, y = standardized_rmse, 
                   group = estimator_name, color = estimator_name)) +
    geom_line(size = 1.2) +
    geom_point(size = 3) +
    geom_hline(yintercept = 1.0, linetype = "dashed", color = "gray30") +
    scale_color_brewer(palette = "Set1") +
    labs(
      title = paste("Standardized RMSE -", ep_params$display_name),
      subtitle = "RMSE relative to 'subgroup (no shrinkage)' estimator",
      y = "Standardized RMSE",
      x = "Simulation Scenario",
      color = "Estimator"
    ) +
    theme_bw() +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      axis.title = element_text(size = 11),
      legend.title = element_text(face = "bold"),
      panel.grid.minor = element_blank()
    )
  
  print(p2)
  cat("\n\n")
  
  # Calculate subgroup performance
  cat("\n### Subgroup-Level Performance Metrics\n\n")
  
  subgroupPerformance <- calculate_subgroup_performance(
    ep_data$merged_df_with_pop, 
    ep_params
  )

  # Focused check: RMSE, bias, coverage for subgroups 4a/4b/4c in scenarios 2 & 3
  focus_keys <- c("4a", "4b", "4c")
  focus_scen <- c("2", "3")
  focus_estimators <- c("subgroup", "Global_r2d2_mid", "Global_r2d2_mid_unshrunk_x4")

  focus_table <- subgroupPerformance %>%
    filter(scenario_no %in% focus_scen,
           join_key %in% focus_keys,
           estimator_name %in% focus_estimators) %>%
    arrange(scenario_no, join_key, estimator_name) %>%
    transmute(
      Scenario = scenario_no,
      Subgroup = join_key,
      Estimator = estimator_name,
      RMSE = round(RMSE, 3),
      Bias = round(bias, 3),
      Coverage = ifelse(is.na(coverage), NA, round(coverage, 3))
    )

  if (nrow(focus_table) > 0) {
    cat('\n#### Focus: Subgroups 4a/4b/4c — Scenarios 2 & 3 (check unshrunk x4)\n\n')
    print(kable(focus_table, caption = paste('Focused metrics -', ep_params$display_name), digits = 3))
    cat('\n\n')
  }
  
  # Generate summary tables
  sum_RMSEbias <- subgroupPerformance %>%
    group_by(estimator_name) %>%
    summarize(
      RMSE = paste(display(mean(RMSE, na.rm = TRUE), 2), " (", 
                  display(min(RMSE, na.rm = TRUE), 2), "-", 
                  display(max(RMSE, na.rm = TRUE), 2), ")", sep = ""),
      abs_bias = paste(display(mean(abs(bias), na.rm = TRUE), 2), " (", 
                      display(min(abs(bias), na.rm = TRUE), 2), "-", 
                      display(max(abs(bias), na.rm = TRUE), 2), ")", sep = "")
    ) %>%
    arrange(estimator_name)
  
  sum_coverage <- subgroupPerformance %>%
    group_by(estimator_name) %>%
    summarize(
      coverage = ifelse(
        all(is.na(coverage)), "NA",
        paste(getPercentage(mean(coverage, na.rm = TRUE)), " (", 
              getPercentage(min(coverage, na.rm = TRUE)), "-", 
              getPercentage(max(coverage, na.rm = TRUE)), ")", sep = "")
      )
    ) %>%
    arrange(estimator_name)
  
  sum_RMSEbias_het <- subgroupPerformance %>%
    group_by(heterogeneous, estimator_name) %>%
    summarize(
      RMSE = paste(display(mean(RMSE, na.rm = TRUE), 2), " (", 
                  display(min(RMSE, na.rm = TRUE), 2), "-", 
                  display(max(RMSE, na.rm = TRUE), 2), ")", sep = ""),
      abs_bias = paste(display(mean(abs(bias), na.rm = TRUE), 2), " (", 
                      display(min(abs(bias), na.rm = TRUE), 2), "-", 
                      display(max(abs(bias), na.rm = TRUE), 2), ")", sep = "")
    ) %>%
    ungroup() %>%
    arrange(estimator_name)
  
  sum_coverage_het <- subgroupPerformance %>%
    group_by(heterogeneous, estimator_name) %>%
    summarize(
      coverage = ifelse(
        all(is.na(coverage)), "NA",
        paste(getPercentage(mean(coverage, na.rm = TRUE)), " (", 
              getPercentage(min(coverage, na.rm = TRUE)), "-", 
              getPercentage(max(coverage, na.rm = TRUE)), ")", sep = "")
      )
    ) %>%
    ungroup() %>%
    arrange(estimator_name)
  
  # Build simplified table - just the summary statistics by estimator
  summary_table <- bind_rows(
    sum_RMSEbias %>% 
      left_join(sum_coverage, by = "estimator_name") %>%
      mutate(
        Estimator = estimator_name,
        RMSE_overall = RMSE,
        Bias_overall = abs_bias,
        Coverage = coverage
      ) %>%
      dplyr::select(Estimator, RMSE_overall, Bias_overall, Coverage)
  )
  
  # Calculate footnote info
  n_scenarios <- n_distinct(subgroupPerformance$scenario_no)
  n_subgroups_per_scenario <- n_distinct(subgroupPerformance$join_key)
  total_subgroups <- n_scenarios * n_subgroups_per_scenario
  
  estimator_names <- sort(unique(subgroupPerformance$estimator_name))
  count_het <- subgroupPerformance %>%
    filter(estimator_name == estimator_names[1]) %>%
    group_by(scenario_no) %>%
    summarise(total_het = sum(heterogeneous))
  total_het <- sum(count_het$total_het)
  total_homo <- total_subgroups - total_het
  
  # Print table
  tbl <- kbl(summary_table,
    caption = paste("Performance Metrics Summary -", ep_params$display_name),
    booktabs = TRUE,
    col.names = c("Estimator", "RMSE (all)", "Abs Bias (all)", "Coverage (95%)")
  ) %>%
    footnote(
      general = c(
        paste("Summary across", total_subgroups, "subgroups (", 
              n_subgroups_per_scenario, "subgroups x", n_scenarios, 
              "scenarios) including", total_homo, "homogeneous and", 
              total_het, "heterogeneous subgroups."),
        "Monte Carlo SE of estimated coverage is ~0.7% given 1000 simulations and true coverage of 95%."
      ),
      general_title = ""
    )
  
  print(tbl)
  cat("\n\n")
  
  # Plot 3: RMSE by Subgroup (Scenario 3)
  cat("\n### RMSE by Subgroup (Scenario 3: Heterogeneous Effects)\n\n")
  
  # Get number of unique estimators for dynamic color palette
  n_est <- n_distinct(subgroupPerformance %>% filter(scenario_no == 3) %>% pull(estimator_name))
  
  p3 <- ggplot(subgroupPerformance %>% filter(scenario_no == 3), 
               aes(x = RMSE, y = fct_rev(join_key), 
                   color = estimator_name, shape = estimator_name)) +
    geom_jitter(size = 4, width = 0, height = 0.2) + 
    scale_colour_brewer(palette = "Spectral", type = "qualitative") +
    scale_shape_manual(values = 1:n_est) +
    labs(
      title = paste("RMSE by Subgroup (Scenario 3) -", ep_params$display_name),
      x = "Root Mean Squared Error",
      y = "Subgroup",
      color = "Estimator",
      shape = "Estimator"
    ) +
    theme_bw() +
    theme(
      plot.title = element_text(face = "bold", size = 14),
      axis.title = element_text(size = 11),
      legend.title = element_text(face = "bold")
    )
  
  print(p3)
  cat("\n\n")
  
  cat("---\n\n")
}

Time-to-Event (Survival) ( AHR )

Overall RMSE by Scenario

Standardized RMSE (Relative to Subgroup Estimator)

Subgroup-Level Performance Metrics

Focus: Subgroups 4a/4b/4c — Scenarios 2 & 3 (check unshrunk x4)

Focused metrics - Time-to-Event (Survival)

Scenario	Subgroup	Estimator	RMSE	Bias	Coverage
2	4a	Global_r2d2_mid	0.301	-0.235	0.770
2	4a	Global_r2d2_mid_unshrunk_x4	0.211	0.006	0.950
2	4a	subgroup	0.203	0.003	0.954
2	4b	Global_r2d2_mid	0.232	0.154	0.893
2	4b	Global_r2d2_mid_unshrunk_x4	0.258	-0.008	0.939
2	4b	subgroup	0.248	0.009	0.944
2	4c	Global_r2d2_mid	0.226	0.148	0.876
2	4c	Global_r2d2_mid_unshrunk_x4	0.253	-0.023	0.949
2	4c	subgroup	0.248	-0.009	0.953
3	4a	Global_r2d2_mid	0.459	0.363	0.704
3	4a	Global_r2d2_mid_unshrunk_x4	0.279	-0.023	0.940
3	4a	subgroup	0.271	-0.004	0.944
3	4b	Global_r2d2_mid	0.209	-0.125	0.909
3	4b	Global_r2d2_mid_unshrunk_x4	0.219	0.005	0.943
3	4b	subgroup	0.211	-0.001	0.945
3	4c	Global_r2d2_mid	0.204	-0.119	0.922
3	4c	Global_r2d2_mid_unshrunk_x4	0.217	0.010	0.950
3	4c	subgroup	0.211	0.003	0.963

Performance Metrics Summary - Time-to-Event (Survival)

Estimator	RMSE (all)	Abs Bias (all)	Coverage (95%)
Global_horseshoe_low	0.16 (0.13-0.49)	0.04 (0.00-0.36)	94% (71%-98%)
Global_horseshoe_mid	0.16 (0.13-0.55)	0.05 (0.00-0.41)	92% (59%-97%)
Global_horseshoe_strong	0.16 (0.13-0.56)	0.05 (0.00-0.43)	92% (56%-97%)
Global_r2d2_low	0.16 (0.13-0.47)	0.04 (0.00-0.38)	95% (66%-99%)
Global_r2d2_mid	0.16 (0.13-0.46)	0.04 (0.00-0.36)	95% (70%-99%)
Global_r2d2_mid_unshrunk_x4	0.17 (0.14-0.39)	0.03 (0.00-0.32)	95% (81%-98%)
Global_r2d2_strong	0.16 (0.13-0.45)	0.04 (0.00-0.36)	95% (74%-99%)
OVAT_full_hierarchical	0.19 (0.14-0.33)	0.02 (0.00-0.18)	96% (89%-99%)
population	0.16 (0.13-0.68)	0.06 (0.00-0.67)	87% (0%-96%)
subgroup	0.23 (0.14-0.36)	0.01 (0.00-0.02)	95% (93%-97%)
Summary across 150 subgroups ( 25 subgroups x 6 scenarios) including 121 homogeneous and 29 heterogeneous subgroups.
Monte Carlo SE of estimated coverage is ~0.7% given 1000 simulations and true coverage of 95%.

RMSE by Subgroup (Scenario 3: Heterogeneous Effects)

---

Binary Outcome ( OR )

Overall RMSE by Scenario

Standardized RMSE (Relative to Subgroup Estimator)

Subgroup-Level Performance Metrics

Focus: Subgroups 4a/4b/4c — Scenarios 2 & 3 (check unshrunk x4)

Focused metrics - Binary Outcome

Scenario	Subgroup	Estimator	RMSE	Bias	Coverage
2	4a	Global_r2d2_mid	0.381	-0.326	0.654
2	4a	Global_r2d2_mid_unshrunk_x4	0.263	-0.007	0.941
2	4a	subgroup	0.259	-0.006	0.952
2	4b	Global_r2d2_mid	0.222	0.144	0.910
2	4b	Global_r2d2_mid_unshrunk_x4	0.273	-0.013	0.944
2	4b	subgroup	0.271	-0.009	0.945
2	4c	Global_r2d2_mid	0.215	0.135	0.905
2	4c	Global_r2d2_mid_unshrunk_x4	0.238	-0.009	0.947
2	4c	subgroup	0.235	-0.007	0.956
3	4a	Global_r2d2_mid	0.480	0.389	0.643
3	4a	Global_r2d2_mid_unshrunk_x4	0.275	0.005	0.937
3	4a	subgroup	0.273	0.008	0.948
3	4b	Global_r2d2_mid	0.258	-0.166	0.850
3	4b	Global_r2d2_mid_unshrunk_x4	0.281	-0.002	0.933
3	4b	subgroup	0.280	-0.003	0.932
3	4c	Global_r2d2_mid	0.244	-0.154	0.840
3	4c	Global_r2d2_mid_unshrunk_x4	0.231	0.000	0.941
3	4c	subgroup	0.228	0.001	0.947

Performance Metrics Summary - Binary Outcome

Estimator	RMSE (all)	Abs Bias (all)	Coverage (95%)
Global_horseshoe_low	0.17 (0.14-0.63)	0.04 (0.00-0.53)	94% (59%-98%)
Global_horseshoe_mid	0.17 (0.14-0.68)	0.05 (0.00-0.60)	92% (45%-97%)
Global_horseshoe_strong	0.17 (0.14-0.69)	0.05 (0.00-0.62)	92% (40%-97%)
Global_r2d2_low	0.17 (0.14-0.59)	0.04 (0.00-0.52)	94% (59%-99%)
Global_r2d2_mid	0.17 (0.14-0.59)	0.04 (0.00-0.52)	94% (61%-99%)
Global_r2d2_mid_unshrunk_x4	0.18 (0.14-0.57)	0.03 (0.00-0.49)	95% (65%-99%)
Global_r2d2_strong	0.17 (0.14-0.61)	0.04 (0.00-0.52)	94% (63%-99%)
OVAT_hierarchical	0.20 (0.15-0.42)	0.02 (0.00-0.26)	96% (91%-99%)
population	0.17 (0.14-0.76)	0.06 (0.00-0.75)	88% (0%-96%)
subgroup	0.25 (0.16-0.39)	0.01 (0.00-0.03)	95% (93%-97%)
Summary across 150 subgroups ( 25 subgroups x 6 scenarios) including 122 homogeneous and 28 heterogeneous subgroups.
Monte Carlo SE of estimated coverage is ~0.7% given 1000 simulations and true coverage of 95%.

RMSE by Subgroup (Scenario 3: Heterogeneous Effects)

---

Count Data ( RR )

Overall RMSE by Scenario

Standardized RMSE (Relative to Subgroup Estimator)

Subgroup-Level Performance Metrics

Focus: Subgroups 4a/4b/4c — Scenarios 2 & 3 (check unshrunk x4)

Focused metrics - Count Data

Scenario	Subgroup	Estimator	RMSE	Bias	Coverage
2	4a	Global_r2d2_mid	0.381	-0.328	0.684
2	4a	Global_r2d2_mid_unshrunk_x4	0.262	0.006	0.959
2	4a	subgroup	0.259	0.000	0.946
2	4b	Global_r2d2_mid	0.235	0.155	0.908
2	4b	Global_r2d2_mid_unshrunk_x4	0.283	-0.017	0.957
2	4b	subgroup	0.275	-0.009	0.944
2	4c	Global_r2d2_mid	0.225	0.149	0.911
2	4c	Global_r2d2_mid_unshrunk_x4	0.252	-0.013	0.955
2	4c	subgroup	0.246	-0.001	0.941
3	4a	Global_r2d2_mid	0.500	0.424	0.657
3	4a	Global_r2d2_mid_unshrunk_x4	0.277	-0.010	0.954
3	4a	subgroup	0.269	0.000	0.949
3	4b	Global_r2d2_mid	0.257	-0.176	0.875
3	4b	Global_r2d2_mid_unshrunk_x4	0.272	0.001	0.942
3	4b	subgroup	0.266	-0.003	0.935
3	4c	Global_r2d2_mid	0.253	-0.174	0.866
3	4c	Global_r2d2_mid_unshrunk_x4	0.239	-0.008	0.947
3	4c	subgroup	0.237	-0.012	0.937

Performance Metrics Summary - Count Data

Estimator	RMSE (all)	Abs Bias (all)	Coverage (95%)
Global_horseshoe_low	0.18 (0.15-0.64)	0.06 (0.00-0.51)	94% (65%-98%)
Global_horseshoe_mid	0.19 (0.15-0.68)	0.07 (0.00-0.58)	92% (52%-97%)
Global_horseshoe_strong	0.19 (0.15-0.69)	0.07 (0.00-0.60)	92% (49%-97%)
Global_r2d2_low	0.18 (0.15-0.59)	0.06 (0.00-0.51)	95% (63%-99%)
Global_r2d2_mid	0.18 (0.15-0.60)	0.06 (0.00-0.51)	95% (66%-99%)
Global_r2d2_mid_unshrunk_x4	0.19 (0.15-0.59)	0.04 (0.00-0.49)	96% (70%-99%)
Global_r2d2_strong	0.18 (0.15-0.62)	0.06 (0.00-0.51)	95% (67%-99%)
OVAT_full_hierarchical	0.21 (0.16-0.47)	0.02 (0.00-0.26)	96% (87%-99%)
OVAT_hierarchical	0.21 (0.16-0.45)	0.02 (0.00-0.23)	96% (87%-99%)
population	0.18 (0.15-0.87)	0.06 (0.00-0.86)	88% (0%-95%)
subgroup	0.26 (0.16-0.44)	0.01 (0.00-0.03)	94% (92%-96%)
Summary across 150 subgroups ( 25 subgroups x 6 scenarios) including 121 homogeneous and 29 heterogeneous subgroups.
Monte Carlo SE of estimated coverage is ~0.7% given 1000 simulations and true coverage of 95%.

RMSE by Subgroup (Scenario 3: Heterogeneous Effects)

---

Continuous Outcome ( Mean Difference )

Overall RMSE by Scenario

Standardized RMSE (Relative to Subgroup Estimator)

Subgroup-Level Performance Metrics

Focus: Subgroups 4a/4b/4c — Scenarios 2 & 3 (check unshrunk x4)

Focused metrics - Continuous Outcome

Scenario	Subgroup	Estimator	RMSE	Bias	Coverage
2	4a	Global_r2d2_mid	0.298	-0.236	0.798
2	4a	Global_r2d2_mid_unshrunk_x4	0.276	0.018	0.947
2	4a	subgroup	0.275	0.018	0.946
2	4b	Global_r2d2_mid	0.204	0.117	0.931
2	4b	Global_r2d2_mid_unshrunk_x4	0.274	-0.002	0.945
2	4b	subgroup	0.273	-0.002	0.949
2	4c	Global_r2d2_mid	0.201	0.110	0.929
2	4c	Global_r2d2_mid_unshrunk_x4	0.244	0.009	0.955
2	4c	subgroup	0.243	0.008	0.953
3	4a	Global_r2d2_mid	0.403	0.343	0.655
3	4a	Global_r2d2_mid_unshrunk_x4	0.277	0.020	0.947
3	4a	subgroup	0.275	0.018	0.946
3	4b	Global_r2d2_mid	0.225	-0.140	0.898
3	4b	Global_r2d2_mid_unshrunk_x4	0.275	-0.004	0.947
3	4b	subgroup	0.273	-0.002	0.949
3	4c	Global_r2d2_mid	0.220	-0.130	0.900
3	4c	Global_r2d2_mid_unshrunk_x4	0.244	0.008	0.953
3	4c	subgroup	0.243	0.008	0.953

Performance Metrics Summary - Continuous Outcome

Estimator	RMSE (all)	Abs Bias (all)	Coverage (95%)
Global_horseshoe_low	0.18 (0.15-0.63)	0.05 (0.00-0.52)	95% (64%-98%)
Global_horseshoe_mid	0.18 (0.15-0.69)	0.05 (0.00-0.59)	93% (52%-97%)
Global_horseshoe_strong	0.18 (0.15-0.70)	0.05 (0.00-0.61)	93% (48%-97%)
Global_r2d2_low	0.18 (0.15-0.59)	0.04 (0.00-0.51)	95% (63%-98%)
Global_r2d2_mid	0.18 (0.15-0.60)	0.04 (0.00-0.51)	95% (64%-98%)
Global_r2d2_mid_unshrunk_x4	0.19 (0.15-0.56)	0.03 (0.00-0.45)	96% (72%-99%)
Global_r2d2_strong	0.18 (0.15-0.62)	0.04 (0.00-0.51)	95% (66%-98%)
OVAT_full_hierarchical	0.20 (0.16-0.25)	0.01 (0.00-0.02)	97% (95%-98%)
OVAT_hierarchical	0.21 (0.16-0.46)	0.03 (0.00-0.29)	96% (87%-99%)
population	0.18 (0.15-0.80)	0.06 (0.00-0.78)	89% (0%-95%)
subgroup	0.26 (0.17-0.41)	0.01 (0.00-0.02)	95% (94%-96%)
Summary across 150 subgroups ( 25 subgroups x 6 scenarios) including 121 homogeneous and 29 heterogeneous subgroups.
Monte Carlo SE of estimated coverage is ~0.7% given 1000 simulations and true coverage of 95%.

RMSE by Subgroup (Scenario 3: Heterogeneous Effects)

---

Summary

This document provides a comprehensive analysis of simulation results across all four endpoint types. Key findings:

1. Cross-Endpoint Consistency: The relative performance of different estimators is generally consistent across endpoint types.

2. Shrinkage Benefit: Global models with shrinkage priors (Horseshoe, R2D2) consistently outperform naive subgroup analysis, especially in scenarios with heterogeneous treatment effects.

3. Scenario-Specific Patterns:

   • Scenarios 1-2 (homogeneous/mild heterogeneity): All methods perform similarly
   • Scenarios 4-5 (random heterogeneity): Shrinkage methods show clear advantages
   • Scenario 6 (misspecified): Highlights importance of model specification

4. Coverage: Bayesian credible intervals maintain appropriate coverage across scenarios and endpoints.

The standardized metrics (log-scale for TTE/binary/count, natural scale for continuous) enable fair comparisons across endpoint types while respecting the natural interpretation of each outcome measure.